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Abstract. In this brief contribution to the Proceedings of the NATO-ASI 
on "Electrostatic Effects in Soft Matter and Biophysics" [1] , which took 
place in Les Houches from Oct. 1-13, 2000, we summarize in short aspects 
of the simulations methods to study charged systems. After describing some 
basics of Monte Carlo and Molecular dynamics techniques, we describe a 
few methods to compute long range interactions in periodic systems. After 
a brief detour to mean-field models, we describe our results obtained for 
flexible polyelectrolytes in good and bad solvents. We follow with a de- 
scription of the inhomogeneity of the counterion distribution around finite 
chains, and continue then with infinitely long, rodlike systems. The last 
part is devoted to the phenomenon of overcharging for colloidal particles 
and its explanation in terms of simple electrostatic arguments. 



1. Introduction 

Polyelectrolytes represent a broad and interesting class of materials [2] that 
enjoy an increasing attention in the scientific community. For example, in 
technical applications polyelectrolytes are wildly used as viscosity modifiers, 
precipitation agents, superabsorbers, or leak protectors. In biochemistry 
and molecular biology they are of interest because virtually all proteins, as 
well as DNA, are polyelectrolytes. 

In contrast to the theory of neutral polymer systems, which is well 
developed, the theory of polyelectrolytes faces several difficulties. Simple 
scaling theories, which have been proven so successfully in neutral polymer 
theory, have to deal with additional length scales set by the long range 
Coulomb interaction [3]. Furthermore there is a delicate interplay between 
the electrostatic interaction of the distribution of the counterions and the 
conformational degrees of freedom, which in turn are governed by a host 
of short range interactions, which renders the problem difficult. There are 
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only two limiting cases that are easy to solve. These are the case of high salt 
excess, effectively screening out the electrostatic interaction (which in turn 
allows one to treat it as a perturbation), or the case of an overwhelming 
dominance of the Coulomb force, which results in a strongly elongated 
chain. Unfortunately it is often just the intermediate case, which proves 
to be the most interesting regime in terms of application, experiment and 
theory. 

Computational simulations provide some unique ways to elucidate the 
properties of charged systems. We first give a more general introduction to 
the relevant simulation methods, and focus then on some recently obtained 
results. 

2. Simulations techniques 

2.1. SOME BASICS ON SIMULATIONS 

There exist a number of nice reviews and books [4, 5, 6] which deal exten- 
sively with various aspects of computer simulations of complex, not neces- 
sarily charged, systems. Thus, the present introduction can also be viewed 
as a guide to the literature. 

There are two basic concepts, which are used in computer simulations 
of complex systems. The conceptionally most direct approach is the molec- 
ular dynamics (MD) method. One numerically solves Newton's equation of 
motion for a collection of particles, which interact via a suitable interaction 
potential U(fi), where are the positions of the particles. Through the 
equation of motion a natural time scale is built in, though, this might not 
be the physically realistic time scale (e. g. if the solvent is replaced by a 
dielectric background). Running such a simulation samples phase space for 
the considered system deterministically. Though this sounds very simple, 
there are many technical and conceptual complications, which one encoun- 
ters on the way. The second approach, the Monte Carlo (MC) method, 
samples phase space stochastically. Monte Carlo is intrinsically stable but 
has no natural time scale built in. This can be reinterpreted, however, by 
an adjustment of suitable "time amplitudes". The MD and MC approaches 
are the basic simulation methods for exploring the statistical properties of 
complex fluids. At present, many applications employ variants thereof, or 
even hybrid methods, where combinations of both are used. Before going 
into detail we ask when is either kind of model appropriate? 

At first sight it is tempting to perform a computer simulation of a 
polyelectrolyte solution where all details of the chemical structure of the 
monomers are included. For instance, the chain diffusion constant D could 
be measured by monitoring the mean square displacement of the monomers 
of the chains. This, however, is tempting only at the very first glance. 
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Even for the fastest computers one would need an exceeding amount of 
computer time. As for all disordered, complex, macromolecular materials, 
(charged) polymers are characterized by a hierarchy of different length and 
time scales, and these length and especially the time scales span an ex- 
tremely wide range [7] . On the atomistic level the properties are dominated 
by the local oscillations of bond angles and lengths. 1 The typical time con- 
stant of about 10 -13 sec results in a simulation time step of 10~ 15 sec. On 
the semi-macroscopic level the behavior is dominated by the overall relax- 
ations of conformation of the objects or even larger units (domains etc). 
The times, depending on chain length, temperature and density, can easily 
reach seconds. To cover that many decades in time within a conventional 
computer simulation is certainly impossible at present. On the other hand, 
it is important to relate the chemical structure of a system to its properties. 



2.2. MOLECULAR DYNAMICS 

MD simulations date back to the early fifties. For a rather complete overview 
about simulations in condensed matter we refer to [4]. Consider a cubic 
box of Volume V = L 3 , containing N identical particles. In order to avoid 
surface effects and (as much as possible) finite size effects, one typically 
uses periodic boundary conditions. The particle number density is given 
by p = N/L 3 . The first simulations employed hard spheres of radius R , 
leading to a volume fraction p v = A/3irR 3 p. Though still used extensively 
for some studies on the glass transition of colloidal systems we focus here 
on soft potentials. 

A key thermodynamic quantity, the temperature, is imposed via the 

equipartition theorem y(rj) = ffcT, m being the particle masses. Note 
that in hard sphere systems temperature only defines a time scale, but is 
otherwise irrelevant. One can find many different soft potentials in the liter- 
ature. However, most widely used is the Lennard-Jones potential U LJ (rij), 
derived originally for interactions of noble gases (Ar, Kr ... ), rij being 
the distance between particle i and j. In its simplest form for two identical 
particles it reads 



U LJ {r l3 ) = 4e 

' ij 



(-) 12 "(-) 6 



(1) 



Usually, a cutoff r c is introduced for the range of the interaction. This 
typically varies between 2.5 a (classical LJ interaction with an attractive 
well of depth e, used for the poor solvent chains in Sec. 4.2) and 2 1 / 6 <r (the 

1 For reactions or to study exited states, the electronic structure is treated explicitely. 
Such methods (Carr-Parrinello simulations, quantum chemistry etc.) are beyond the 
scope of the present paper. 
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potential is cut off at the minimum, leading to a pure repulsive interaction, 
as is typically done for good solvent chains). For chain molecules a bonding 
interaction for r < Rq 

1 r 2 
CfeneM = --kRl ln(l - (2) 
I R Q 

is added, which keeps the bond length below a maximum of Rq. The spring 
constant k varies between 5 and 30 e/a 2 . Electrostatics is included via the 
Coulomb interaction 

Ucirij) = — q iqj (3) 

' ij 

e 2 

Here Ib ■= 47r£o£ " fcaT , where eo is the elementary unit charge, ks is the 
Boltzmann constant, T denotes temperature, £o and e r are the vacuum 
and relative dielectric permeability of the solvent, respectively, and qij are 
charges measured in units of eo- Two monovalent charges separated by the 
Bjerrum length Ib have an interaction energy equal to ksT. The Bjerrum 
length thus is a measure of the interaction strength. It is equal to 7.14 A for 
water at room temperature. The computational aspects of the long range 
potential will be discussed shortly in Sec. 2. 4. 

The unit of energy is e, of length a and of mass m. This defi nes the "LJ- 
units" for temperature [T] = e/(kBT), time [t] = \J a 2 m/e and number 
density [p] = a~ 3 . In most practical programs a, m, e are used as the 
basic units and set to one. The straight forward simulation technique is to 
integrate Newton's equations of motion for the particles: 

mik = -V U ( r ij) ( 4 ) 

Since energy in such a simulation is conserved we have a microcanon- 
ical ensemble. Presently other thermodynamic ensembles are commonly 
used for practical applications (NPT: isobaric-isothermal, NVT: isother- 
mal (canonic) ... ). Because we often employ a stochastic MD method, 
known as the Langevin thermostat [8], we will briefly describe its main in- 
gredients. Instead of integrating Newton's equations of motion, one solves 
a set of Langevin equations 

mik = -V U ( r ij) ~ r ^ + ( 5 ) 

with £i(t) being a 5-correlated Gaussian noise source with its first and 
second moments given by 

<£(*)> =0 and {£ l {t)-i J {t')) = Qk B TT8 l3 5{t-t'). (6) 
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The friction term — Tfi and the noise £i(t) is thought of as imitating the 
presence of a surrounding viscous medium responsible for a drag force and 
random collisions, respectively. The second moment of is adjusted 
via an Einstein relation in order to reach the canonical state in the limit 
t — > oo. The dynamics generated by the Langevin equation can alternatively 
be written as a general Fokker-Planck process. This permits a transparent 
proof of two important facts: (i) the stationary state of the process is the 
Boltzmann distribution and (ii) the system will evolve and converge to 
the Boltzmann distribution [9]. Since for small times the stochastic part 
is more important than the deterministic one [y/t 3> t for small t) it is 
actually not necessary to use Gaussian random variables in the simulation 
[10]. It suffices to use equidistributed random variables with first and second 
moment being identical to the Gaussian deviate. 

A simple but very efficient and stable integration scheme is the Verlet 
algorithm (more complicated methods, which do not have time inversion 
symmetry, do not, in general, perform significantly better). With a simula- 
tion time step St , where St « 2-K/uj max and u> ma x is the typical highest 
frequency of the system (for crystals the Einstein frequency), we have in 
one dimension 

5t 2 St 3 

n{t + 5t) = n(t) + Stvi(t) + — Oi(t) + — di(t) + 0{St A ) 

St 2 St 3 

nit- St) = r i (t)-5tv i (t) + —a i (t)-—a i (t)+0(5t 4 ), 
where vi(t) = fi(t) and dj(i) = Vi(t). An addition of the two lines yields 

n{t + St) = 2n(t) - n(t - St) + st 2 ai (t) + o{St 4 ) (7) 

Therefore the position calculations have an algorithmic error of 0(St 4 ). 
Subtraction of the lines yields 

vi(t) = ^[n(t + St) - n (t - St)} + 0(St 3 ) (8) 

leading to errors of 0(St 3 ). There are many variants of this basic integrating 
scheme used throughout the literature [6] . One can follow the realistic time 
evolution of a system, as long as the forces/potentials are realistic for the 
modeled system and as long as classical mechanics is sufficient. In a purely 
deterministic simulation the accumulation of small errors can cause signif- 
icant deviations from the real trajectory. If the system is ergodic, which 
requires mixing of normal modes (recall the well-known Fermi-Pasta-Ulam 
problem, where one asks how anharmonic a potential has to be in order 



6 



to equilibrate a one dimensional chain of particles [11]) one can determine 
ensemble averages from time averages 



of any physical quantity A of interest. This describes the most elementary 
Ansatz for a microcanonical simulation [6]. Here all extensive thermody- 
namic variables of the system, namely N, V, E are kept constant. Some- 
times this is also called NVE ensemble. As mentioned before, most appli- 
cations employ other ensembles such as the canonic (NVT), the isobaric- 
isothermal (NPT) or even the grand canonical (/U, P, T) ensemble, being 
the chemical potential. As a general rule, in cases such as two phase coexis- 
tence or calculations of transport properties, one should choose an ensemble 
with many intensive variables kept constant as possible. For charged sys- 
tems, however, it is rather difficult to perform efficient simulations in the 
(TV, P, T) or (/x, P, T) ensembles. Therefore the most common ensemble is 
the NVT since it is easy to use the deterministic equations with additional 
stochastic terms to constrain the temperature (Eq. 5). 

2.3. MONTE CARLO METHOD 

The classical Monte Carlo approach goes to the other extreme, namely to 
purely stochastic sampling. Starting from a particular configuration, ran- 
domly a particle (or a number of particles) is selected and displaced by a 
random jump. For hard sphere systems the move is accepted if the new con- 
figuration complies with the excluded volume; if not the old configuration 
is retained. The approach is also called simple sampling. This cycle is re- 
peated over and over. Once every particle on average has a chance to move, 
one Monte Carlo step is completed. This is the most basic Monte Carlo 
simulation (see e.g. [4, 12]). Since there is no energy involved, it trivially 
fulfills detailed balance 



where PF({x} — > {y}) is the probability to jump from state {x} to state 
{y} and P eq ({x}) the equilibrium probability of state {x}. All states have 
exactly the same probability. Detailed balance is a sufficient condition for 
a MC simulation to relax into thermal equilibrium, though this may take 
a very long time. Special cases of algorithms without detailed balance will 
not be discussed here. 

It is useful to compare the basic aspects of MC simulation to the exam- 
ples discussed above for molecular dynamics simulations. The Hamiltonian 




(9) 



i=l 



W({x} -+ {y})P eq ({x}) = W({y} -+ {x})P eq ({y}) 



(10) 
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depends for simplicity only on the positions of all particles {r^} and is 
denoted by H({ri}. The expectation value of any observable A is given by 

<A>=Y,A{{n})P eq {{ n }) (li) 

{n} 

with 

Pe q {{n}) = eM-H({ri})/k B T)/Z 

Z = Y, ex P(-H/k B T) (12) 

{n} 

An exact way would be to sample all possible states, which in all but the 
most trivial cases is impossible. Thus we sample phase space stochastically. 
Taking a particle at random, calculating its energy, one moves it and cal- 
culates the new energy. With P({x}) being the Boltzmann-probability of 
the original state and P({y}) being that of the new state, detailed balance 
is obeyed if 

^| 2 {I}) = eM-(H({x}) - H({y})/k B T)} (13) 

Under this condition the algorithm is ergodic and the system relaxes into 
equilibrium. The Metropolis method is the most frequently used prescrip- 
tion one to accept or reject a move: 

w{{x} ^ {y}) = r j exp(tf( W) -H( { y})/k B T) , AH > (M) 

Since only the ratio of the W's is relevant, T is an arbitrary constant be- 
tween zero and one, usually T = 1. A random number x, equally distributed 
between and 1, is used to decide upon the acceptance of a move. If 
x < W({x} — > {y}) the move is accepted, otherwise rejected. (For T = 
1 any move, which lowers the energy is accepted.) This is the basic MC 
procedure used for sampling phase space in statistical physics. 

In many cases, however, one also would like to gain information on 
the dynamics of a system or even better, of a model system. How can a 
MC simulation, with no intrinsic time scale, be used to obtain informa- 
tion on the dynamics? In the method described above the system evolves 
from one state to another by a local move. Through these local stochastic 
moves the configurations of particles change with " time" . This is a dynamic 
MC method based on a Markov process, where subsequent configurations 
{x} — ► {x 1 } — > {x n } — > . . . are generated with a transition probability 



8 



VF({x 1 } — > To a large extent the choice of the move is arbitrary, as 

long as one can interpret it as a local elementary unit of motion. The prefac- 
tor r can actually be interpreted as an attempt rate T = r" 1 for the moves 
and introduces a timescale. This "changes" the purely statistical transition 
probability into a transition probability per unit time[5, 12]. To compare 
simulated (overdamped) dynamics with an experiment, it essentially re- 
quires determination of r G (e. g. diffusion constants). It is obvious, that 
this simulation does not include any hydrodynamic effects since there is no 
momentum involved. There are very interesting, more advanced methods 
like DPD (dissipative particle dynamics) and Lattice-Boltzmann methods, 
currently under development in order to include this efficiently [13]). Using 
the interpretation of a MC step as a time step, ensemble averages can be 
written as time averages: 

1 M If* 

< A >= — — Y A({x 1 }) / dt'A(t') . (15) 

M - M ^ t-t J t v ' K ' 

° i=M a +l ° Jt ° 

We view one attempted move per system particle as one time-step. 

The first configurations in a simulation are usually not yet equilibrium 
configurations. One first has to "relax" the system into equilibrium, mean- 
ing the data for the first M Q steps are omitted. In this interpretation the 
dynamic Monte Carlo procedure is nothing but a numerical realization of 
a Markov process described by a Markovian master equation 

j t P({x},t) = -Y,W{{x} ^ &})P{ix},t) 

+ J2w({x i }^{x})P({x i },t) (16) 

{x*} 

with P({x},t) the time dependent probability of state {x}. The condition of 
detailed balance is sufficient that P eq ({x}) is the steady-state solution of the 
master equation. If all states are mutually accessible P({x},t) must relax 
towards P eq ({x}) as t — > oo irrespective of the starting state. Note however 
that the choice of a "good" starting state can save enormous amounts of 
CPU time. 

So far, the two extreme cases for classical, particle based computer sim- 
ulations were discussed, microcanonical MD and canonical MC. There are 
many approaches "in between" which are used depending on the problem 
under consideration. The techniques range from pure MD, where Newton's 
equations of motion are solved (x = —VU), MD coupled to a heat bath and 
added friction ("Langevin MD", "Noisy MD"), (x = -W - (x + /(*)), C 
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friction, f (t) random force), Brownian Dynamics (BD) (x = —VU + ran- 
dom displacement), force biased MC (attempted moves are selected from 
the very beginning according to local forces), to plain MC as described 
above. 

For the application to polymers one should keep in mind, that the con- 
formational entropy of the chains add additional complications, which make 
proper equilibration especially difficult or time consuming. Thus, wherever 
possible, methods should be used, which are faster than the slow intrinsic 
dynamics of the chains [5]. 

2.4. METHODS FOR LONG RANGE INTERACTIONS 

One of the biggest problems for the simulations of charged systems is the 
long range nature of the Coulomb interactions. In principle, each charge 
interacts which all others, leading to a computational effort of 0(N 2 ) al- 
ready within the central simulation box. For many physical investigations 
one wants to simulate bulk properties and therefore introduces periodic 
boundary conditions to avoid boundary effects. The standard method to 
compute the merely conditionally convergent Coulomb sum 



where the prime denotes that for n = the term i = j has to be omitted, is 
the traditional Ewald summation [14] . The basic idea is to split the original 
sum via a simple transformation into two exponentially convergent parts, 
where the first one, fa, is short ranged and evaluated in real space, the 
other one, fa, is long ranged and can be analytically Fourier transformed 
and evaluated in Fourier space: 



dt, though other choices are possible and sometimes more advantageous[15, 
16, 17]. For any choice of the Ewald parameter a and no truncation in the 
sums the formula yields the exact result. In practice one wants to cut off 
the infinite sum at some finite values r c and k c to obtain E to a user 
controlled accuracy, which is possible by using error estimates [18]. The 
aformentionend procedure results in the well known Ewald formula for the 
energy of the box 




(17) 




E = E^ + E^ + E^ + E^ 



(19) 
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where the contributions from left to right are the real space, Fourier space, 
self , and dipole-correction energy terms. These are given by 

^•>^EE«« (») 
E(k) = ^Er"" /4a2 i^' 2 ( 21 ) 

EiS) = "^E^ 2 ( 22 ) 

* i 

Eli > = orlw(?" r ') 2 ' <23) 

and the Fourier transformed charge density p(k) is defined as 

p(k) = / d 3 r /9(r)e- ik ' r = Y^qj e"^. (24) 

The dipole term depends not on a, hence is independent of the splitting 
function. It reflects the way Eq. (17) is summed up, here in a spherical way 
towards infinity [19]. It also includes a correction for the dielectric constant e' 
outside the summed up sphere volume. For metallic boundary conditions, 
e' = oo, and the dipole term vanishes. In principle the thermodynamic 
properties should be independent of the choice of boundary conditions [20]. 
The Ewald sum has complexity 0(N 3 / 2 ) in its optimal implementation [21], 
and therefore is not suitable for the study of large systems (N > 0(1000)). 
Implementing a fast Fourier transformation (FFT) for the Fourier part 
results in the so-called particle-mesh-Ewald formulations, which improve 
the efficiency to 0(N log N) [22, 23, 24, 25] , and which can also be efficiently 
be parallelized[26, 27]. The most versatile and accurate method of all mesh- 
methods is the oldest P3M algorithm [22, 28], for which also precise error 
estimates exist [29]. Another way of computing Eq.(17) is via a convergence 
factor 

E=]im-TT m > eXP { ~ m ' +HLl) . (25) 
^02^^ \ra +nL\ y ' 

n ij J 

This approach is used in the Lekner [30] and Sperb [31] methods to effi- 
ciently sum up the 3D Coulomb sum. Although the method in its original 
versions has 0(N 2 ) complexity, Sperb et al. have developed a factorization 
approach which yields an 0(N log N) algorithm [32]. 
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Other advanced methods of 0(N log N) are tree algorithms [33], which 
are the first order approximation of even better, so-called fast multipole 
methods [34]. These can reach a linear complexity, but at the expense of 
a heavy computational overhead which makes these methods advantageous 
only for a very large number of charges (N ^ 100 000) [35]. 

For thin polyelectrolytes films or membrane interactions one is also in- 
terested in summations where only 2 dimensions are periodically replicated 
and the third one is of finite thickness h (2D + h geometry). For this geome- 
try Ewald based formulas are only slowly convergent, have mostly 0(N 2 ) 
scalings and no "a priori" error estimates exist [36] . Recently Arnold [37, 38] 
developed a method which is based on convergent factors, whose errors 
are well controlled, and which uses a factorization approach resulting in 
an 0(N 5 / 3 ) scaling (MMM2D). In two dimensions the convergence factor 
based methods and the Ewald sum methods yield exactly the same results, 
there is no dipolar correction term needed [37]. This is in contrast to the 3D 
methods[19]. However, an even better scaling can be achieved, if one returns 
to the 3D Ewald formula, and allows for a large empty space between the 
unwanted replicas in the third dimension [39]. However, so far the method 
has been only checked on a trial and error basis. We recently improved this 
situation by computing an analytic error term which accounts for the con- 
tributions of the unwanted replicas. By simply subtracting this term from 
the 3D sum, which is a linear operation in TV, one can in principle come as 
close as desired to the real 2D + h sum, by allowing for just some arbitrary 
small amount of empty space between the layers [40]. Using then again the 
P3M method we obtain an N log N scaling with well controlled errors also 
for the 2D + h geometry, which up to now seems to be the optimal choice. 

To simulate the structure of water (or other dipolar solvents), one also 
needs to treat the dipolar interactions in a similar fashion. Also here the 
Ewald method is applicable [21], and error estimates exist [41]. 

3. Mean-field models: Debye-Hiickel chains 

While we deal here mainly with systems where the charges are explicitly 
taken into account, historically and even up to now many studies consider 
the ions solely in a mean field approximation. In the first step all non 
bonded charges are considered as a smeared continuous charge density. 
Such a situation is described by the Poisson-Boltzmann (PB) equation. 
While this in most cases is not exactly solvable, many studies employ the 
Debye-Hiickel approximation, which is the solution of the linearized PB 
equation[42, 43, 44, 45]. The resulting potential between charges is the 
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screened Coulomb potential 

V DH = ^T^^, (26) 
e r 

withl/K being the Debye screening length. The polymer can now easily 
be modeled as a random walk of N monomers from which a fraction / is 
monovalently charged. If one in addition introduces the stiffness along the 
backbone of the chain by a cosine-potential, the total Hamiltonian reads 

N-l N i-1 , s 

JL _ . A £ (5 . + gg %( , j)iB ?to). (27) 

The positive amplitude A defines the strength of the angular potential 
and the Heavyside ^-function is one if the monomers i and j are charged, 
and zero otherwise. The symbol 6j is the bond vector between monomer 
number i and i + For the present simulations the bond length |6| and the 
Bjerrum length lg are fixed to one a, and thus set the basic length scale. 
One bond mimics several neutral monomers as usual for coarse grained 
simulation models. The longest chains considered contained up to 2049 
repeat units with a charge fraction of / = jq. Mapping this onto a PS- 
NaPSS copolymer for a more flexible case and keeping in mind that for 
these polymers roughly three repeat units fit into one Bjerrum length one 
arrives at a molecular weight of more than 600 000 g/mol. Thus, for these 
kind of questions computer simulations are quite capable of covering the 
experimentally interesting regime. This allows us to systematically vary not 
only the chain length and the screening length but also the chain bending 
stiffness. This is only of limited experimental relevance since very long 
isolated chains in dilute solutions cannot be experimentally analyzed so 
far. However one of the central questions in the theory of polyelectrolytes 
is whether the characteristic electrostatic length is a linear or quadratic 
function of the Debye length (i.e. proportional to the square of charge 
density or to the charge density itself) . Analytic results mainly predict a k 2 
asymptotic dependency of the persistence length [3]. However, an estimate 
of the required the chain length to reach the asymptotic regime for flexible 
weakly charged polyelectrolytes shows that the chains are so long that 
all experimentally realizable concentrations are in the semi-dilute regime. 
Thus, this is a typical question of interest which is of no direct importance 
for experiments. Nevertheless, it can add significantly to our understanding 
which in turn can also influence the interpretation of experiments in the 
semi-dilute regime. Therefore, it is worthwhile to undertake some efforts 
in computer simulation to investigate this question. There are a number 
of attempts to do this, which so far do not lead to a clear-cut answer. A 
typical result from such a simulation is given in Fig. 1. 



13 




Figure 1. Electrostatic persistence length as a function of l//t for several A values. 



It shows the dependence of the electrostatic persistence length on the 
screening length. The two dotted lines indicate the variational (Var) and the 
Odijk-Skolnik-Fixman (OSF) results ( L e oc ft -1 , or k~ 2 respectively). As 
one can see for the regime covered there seems to be a continuous transition 
from a very weak dependency on k towards the asymptotically expected 
k~ 2 behavior for the semi-flexible polyelectrolyte. In the limit of very stiff 
chains the OSF result becomes exact and the data has to agree with that. 
Even though the simulated chains lie very well within the experimentally 
relevant regime, there is, at least for typical experimental flexibilities, no 
clear sign of the crossover into the asymptotic OSF regime. This allows 
for two important conclusions. First of all, the chains are certainly not 
long enough to display asymptotic behavior. Secondly, typical experimen- 
tal chains are also not long enough to display the asymptotic behavior as 
predicted by mean- field theories. Unlike neutral polymers, polyelectrolytes 
normally are not in an asymptotic limit where the predicted scaling laws 
can be cleanly observed. Another consequence is revealed in a closer analy- 
sis of the simulation data. For semiflexible polyelectrolyte chains there is no 
unique persistence length anymore, as all theoretical pictures assume. Over 
short distances the intrinsic stiffness dominates and gives a clear signal in 
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the bond direction correlation function. Only over large distances along the 
backbone of the chain does the electrostatic effect show up and introduces a 
second characteristic length scale into the system. This can lead to the spe- 
cial situation that the electrostatic contribution to the persistence length 
for stiff chains is smaller than for flexible chains. The reason simply lies in 
the fact that for the stiff chains the charges along the backbone are already 
much further apart because of the intrinsic stiffness as opposed to flexi- 
ble chains. These larger distances then experience a dramatically weaker 
electrostatic repulsion due to the screening of the electrostatic interaction 
resulting in a weaker effect on the persistence length. For details we refer 
to Ref. [46]. 

Altogether, the main conclusion from these simulations of simplified 
mean-field like models is on the one hand that polyelectrolytes have to be 
extremely long to be in the asymptotic regime, and that one has to be 
very careful in deriving general statements from either simulations, which 
typically are not asymptotic just as experiments, and analytic theory, which 
is usually in an asymptotic limit. 

Another conclusion is that the concept of persistence length for poly- 
electrolytes certainly is not well defined in terms of the properties of the 
chains. As soon as intrinsic stiffness is included there is no longer a unique 
length scale that describes the internal structure of the chain. On the other 
hand this is the basis of all classical models used in analytic theory. 

4. Solutions of flexible polyelectrolytes 

4.1. GOOD SOLVENT CHAINS WITH EXPLICIT COUNTERIONS 

The first investigation of totally flexible many chain polyelectrolyte systems 
in good solvent with explicit monovalent counterions was performed several 
years ago [47]. The simulations were carried out mostly with systems of 8 
or 16 chains with N m = 16, 32 and 64. Instead of the P3M algorithm 
a spherical approximation in a truncated octahedral simulation box was 
used which, for values smaller than N tota i « 500 is faster than the PME 
method. More details of the whole study can be found in [47]. All beads and 
counterions interacted with the truncated Lennard-Jones potential plus the 
full Coulomb interaction. 

In this work experimental values of the osmotic pressure and the max- 
imum position in the interchain structure facture were successfully repro- 
duced. One of the important findings was that the chains essentially are 
never rodlike. Counterion-chain correlations can dramatically shrink the 
polyelectrolyte chain. The end-to-end distance shortens significantly as the 
density increases from dilute towards the overlap density. The chain struc- 
ture is highly anisotropic in the very dilute limit, and the scaling with 
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respect to N m is asymmetric; but as the overlap density is approached, 
the structural anisotropy dissipates and the scaling becomes approximately 
symmetric. On long length scales the chain structure continuously changes 
from very elongated to neutral-like coils. Yet, on short length scales, the 
chain structure is density independent and elongated more than neutral 
chains. 

It was found that in the dilute limit the scaling for the extension per- 
pendicular to the chain was R± oc j\r ' 65-0 ' 70 , and for the extension parallel 
R\\ oc jV°' 90 ~ L0 °. Near the density, where the rodlike chains in disordered 
solution overlap, p ~ iV -2 , R± grows at the expense of i?y until at the over- 
lap density p* the effective exponent is about 0.82. The transition regime 
ranges from p ~ N~ 2 to about p ~ N~ 1A where the coils start to overlap 
and one eventually reaches v = 1/2 in the semidilute regime. The expo- 
nents reported should not necessarily be taken as asymptotic (N — ► oo), 
however they should be relevant for many experimental systems. 

4.2. POOR SOLVENT CHAINS WITH EXPLICIT COUNTEPJONS 

Many polyelectrolytes possess a carbon based backbone for which water 
is a poor solvent. Therefore, in aqueous solution, there is a competition 
between the solvent quality, the Coulombic repulsion, and the entropic de- 
grees of freedom. The conformation in these systems can under certain 
conditions assume pearl- necklace like structures [48]. These also exist for 
strongly charged polyelectrolytes at finite densities in the presence of coun- 
terions[49, 27]. The simulations in Ref.[49] used 16 chains of length N m = 
94, with a charge fraction of / = 1/3, and monovalent counterions. The hy- 
drophobic interaction strength was tuned by means of the Lennard-Jones 
parameter e. There we showed that the polymer density p can be used as 

a very simple parameter to separate different conformation regimes. This 

R 2 

can already be seen in the plots of the end-to-end distance R e and r = 

r g 

versus p in Fig. 2. At very high densities the electrostatic interaction is 
highly screened, so that the hydrophobic interaction wins, and the chains 
collapse to dense globules. If one slightly decreases the density, the chains 
can even contract further, because there are no more steric hinderences 
from the other chains or counterions, and the screening is smaller. The 
collapsed globules, however, have still a net charge, and repel each other, 
so that this phase resembles a charged stabilized colloid or microgel phase. 
With decreasing density the electrostatic interaction will dominate over the 
hydrophobic one. The chains will tend to elongate, assuming pearl-necklace 
conformations, like in Fig. 3, as they have been predicted for weakly charged 
polyelectrolytes in Ref. [48]. The more the chain stretches, the smaller be- 
come the locally compact regions. Note that in contrast to the analytical 
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Figure 2. Re (left) and r (right) versus density p for hydrophilic (phil), weak hydrophobic 
(e_tj = 0.5), and strongly hydrophobic (etj = 1.5) chains 




Figure 3. Typical polyelectrolyte conformation for a density p = 2 • 10 4 a 3 , showing 
5 pearls. The chain had 382 monomers with a charge fraction / = 1/3, Ib = 1.5, and 
e=1.75. 

theories[50, 51], the pearls are stable, even though there are counterions 
localized near and/or inside the pearls. 

Experimentally there are some hints for the existence of pearl-necklace 
chains[52, 53]. One of the obstacles to observing them in scattering ex- 
periments could be related to the strong fluctuations of the pearl number. 
Even in equilibrium we have found coexistence of several pearl states[27]. 
In Fig. 4.2 we see the time evolution of one single chain composed of 382 
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Figure 4- Time development of the radius of gyration Rq, the hydrodynamic radius 
Rh and the number of observed pearls for the same system as in Fig. 3. 



monomers with a charge fraction / = 1/3, Ib = 1-5, and e = 1.75 in a many 
chain system at density p = 1.48 x 10 -5 . One observes jumps between a 
five and four pearl configuration. Also the position of the pearls move quite 
vividly 

The different length scales appearing in a chain can be analyzed by 
looking at the spherically averaged form factor S\(q) of the chain. The 
maximum seen at q = 6 comes from the monomer extension. In the range 
1 < q < 2 we observe a sharp decrease in Si, which comes from the scat- 
tering from the pearls, because it shows the typical Porod scattering of 
Si(q) — q~ 4 - The kink at q « 1.66 appears at the position expected from 
the pearl size, but is broadly smeared out due to large size fluctuations. The 
shoulder which can be seen at q ~ 0.5 does not come from the intra-pearl 
scattering but is due to the scattering of neighboring pearls along the chain 
(inter-pearl contribution), which have a mean distance of (rpp) = 13. 3. a. 
It is also smeared out due to large distribution of inter-pearl distances. 
We conclude that the signatures of the pearl-necklaces are weak already 
for monodisperse samples. A possible improvement could be achieved for 
chains of very large molecular weights and only few pearl numbers, which 




could lead to stable and large signatures. Many more interesting results on 
poor solvent polyelectrolytes can be found in Ref.[27] and will be published 
soon. 

4.3. COUNTERION DISTRIBUTION AROUND FINITE 
POLYELECTROLYTES 

We recently completed a study of the spatial distribution of the counterions 
around strongly charged, flexible polyelectrolytes [54] in good and poor sol- 
vent. There we demonstrated that by partially neutralizing the quenched 
charged distribution on the chain backbone the inhomogeneous distribu- 
tion of counterions leads to the same qualitative effects that are observed 
in weakly charged polyelectrolytes with an annealed charge distribution [3]. 
This is due to the presence of the mobile partially neutralizing counterions, 
which results in an annealed backbone charge distribution. The common 
underlying physical mechanism for the end-effect is the differences in the 
electrostatic field of the chain along its backbone. The strength of the end- 
effect depends on parameters like chain length, charge fraction and ionic 
strength, and those dependencies were found in agreement with the scaling 
predictions. We found a saturation of the end-effect for long chains, when 
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the chain extension, namely R e , is at least twice as large as the Debye 
screening length. A simple Debye-length criterion appeared to be sufficient 
to explain the penetration depth of the end-effect. However, looking at the 
amplitude dependency on density and ionic strength of the solution, we 
found that both parameters, the number of annealing ions and the ionic 
strength of the solution, influence the end-effect and that the first one dom- 
inated. The amplitude of the end-effect was shown to depend strongly on 
the charge parameter £ := Is/b, where b is the distance of the bare charges 
on the backbone of the chain. The definition of such an end-effect via close 
mobile counterions can not be made for an effective charge £ << 1, because 
under dilute conditions there are almost no counterions close to the chain. 

Even though the chain conformation is very different in the poor solvent 
case the end-effect was found to be qualitatively the same, namely the 
counterions are more likely to be found at the middle of the chain than 
at the ends. We could also clearly see the necklace structure by looking at 
the effective charge along the contour length. However, the string length of 
our simulated pearl-necklaces was too short to show any charge difference 
between the pearls and the strings, as has been predicted in Ref. [55]. 

We also obtained a fairly good agreement of the simulated ion distribu- 
tion with the PB solution of the cell model of an infinitely extended charged 
rod [42]. This supports the idea that the description of polyelectrolytes as 
rodlike objects in mean-field theory is valid in the dilute regime. Further 
improvements could probably be achieved along the lines of Ref. [56], where 
a combination of a cylindrical and spherical cell model is used to describe 
the solution properties of polyelectrolytes. 

4.4. RODLIKE POLYELECTROLYTES 

Stiff linear polyelectrolytes can be approximated by charged cylinders. 
This is a relevant special case, applying to quite a few biologically im- 
portant polyelectrolytes with a large persistence length, like DNA, actin 
filaments or microtubules. Within PB theory [57] and on the level of a cell 
model the cylindrical geometry can be treated exactly in the salt-free case 
[42, 58, 59, 60, 61, 62, 63], providing for instance new insights into the phe- 
nomenon of the Manning condensation [64, 65]. For low line charges, the 
agreement between PB theory and the simulations of the full interaction 
system is rather nice. However, PB theory fails quantitatively (underes- 
timated condensation) and qualitatively (overcharging, charge oscillations 
and attractive interactions); see, e.g. Ref. [63, 66, 67]. 

Recently the osmotic coefficient of a synthetic stiff polyelectrolyte, a 
poly(para-phenylene), was measured in a salt-free environment [68, 69]. We 
have compared this data to predictions of PB theory, and a local density 
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functional theory which includes a correlation correction of the basis of a 
recently proposed Debye-Hiickel-Hole-cavity theory (DHHC) [70], and sim- 
ulational results within the cell model. We find that correlation effects en- 
hance condensation and lower the osmotic pressure, yet are not fully able to 
explain the discrepancy with the experimental data. Here the approach of 
working within the "primitive model" breaks down. In our opinion, specific 
interactions between the counterions, the macroion, and the solvent par- 
ticles are needed to explain the discrepancy. Other theoretical approaches 
beyond the cell model which try to incorporate finite-size effects and in- 
teractions of the macroion itself will in general lead to a higher osmotic 
coefficient which is in contrast to the experimental data[71]. 

Attractive interactions have been observed[72, 73, 74, 75] and predicted 
between like-charged macromolecules. However, there are nice rigorous re- 
sults which prove that these effects cannot be described by mean-field 
theories [76, 77, 78]. Especially in the community of biological inspired 
physics[45, 79, 80, 81], these interactions are thought to be important for 
the clarification of the mechanism behind DNA compactification in viral 
heads [82], the chromatin structure [83], and novel methods for gene deliv- 
ery [84], to name just the most prominent examples. There are numerous 
simulations which show similar attractions on a distance of few counterion 
diameters [66, 85, 86, 87, 88, 89, 90]. 

The mechanism which is driving the observed attractions for rod-like 
systems has been speculated to be correlations between the counterion lay- 
ers around the macroion. However, until now, no unique theoretical picture 
has emerged that can clarify the detailed mechanism behind the attrac- 
tions. There is the low temperature Wigner crystal theory, initiated by Refs. 
[91, 92, 81], which postulates an ordered ground state of the counterions. 
Then there are theories which are based on Van der Waals type correlated 
fluctuations [60, 93, 94, 95], that are in principle hight T theories. There are 
also theories which are fluctuation based, but are valid at low T[44, 96, 97]. 
Integral equation [43, 98, 99] theories on various approximation levels have 
been demonstrating the existence of these attractions for a long time, but 
from these theories it is difficult to extract the detailed mechanism behind 
the observed correlations. Here also simulations can be helpful, because they 
have in principal access to all correlations [100]. More details of our results 
in rod-like geometries can be found in Refs. [63, 66, 67, 70, 71, 101, 102]. 

5. The energetic path to understand overcharging 

There has been a recent interest in the study of systems which are strongly 
coupled by Coulomb interactions. These systems show a variety of, at first 
sight, surprising behaviors, which can not be accounted for by the mean- 
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field PB theory. For example, there are attractions between like charged 
objects and a charge reversal of macroions occurs when viewed from some 
distance. This means that there are more ions of the opposite charge within 
a certain radius around the macroion then necessary to charge neutralize 
it. This overcompensation is called "overcharging". 

In this section we want to demonstrate that there are situations for 
charged colloidal objects in which one can understand the phenomenon 
of overcharging by very simple energetic arguments. By overcharging, in 
general, we mean that the bare charge of the macroion is overcompensated 
at some distance by oppositely charged "microions". To achieve this in 
nature we have to add salt to the system. For the sake of simplicity, however, 
we will consider non-neutral systems, because they can on a very simple 
basis explain why colloids prefer to be overcharged. 

5.1. THE MODEL 

Our model is solely based on electrostatic energy considerations, meaning 
that we only look at the ground state of a system of charges. We con- 
sider a colloid of radius a with a central charge Z. In the ground state the 
counterions of this colloid are located on the surface, because there they 
are closest to the central charge. On the other hand they want to be in 
such a configuration that they minimize their mutual repulsion. For two, 
three, and four counterions these configurations correspond to a line, an 
equilateral triangle, and a tetrahedron, respectively, regardless of the cen- 
tral charge magnitude. The problem of the minimal energy configuration 
of electrons disposed on the surface of a sphere dates back to Thomson 
[103], and is actually unsolved for large TV. The reason is, that there are 
many metastable states which differ only minimally in energy, and their 
number seems to grow exponentially with N. Also chemists developed the 
valence-shell electron-pair repulsion (VSEPR) theory [104] which uses sim- 
ilar arguments to predict the molecular geometry in covalent compounds, 
also known as the Gillespie rule. 

A simple illustration of energetically driven overcharging is depicted in 
Fig. 6. The central charge is +2e, and the neutral system has two coun- 
terions of valence 1. If we add successively more counterions of the same 
valence, and put them on the surface such that their mutual repulsion is 
minimized, we can compute the total electrostatic energy according to 

E(n) = k B T{l B /a) [-nZ m + , (28) 

where /(#?) is the repulsive energy part which is only a function of the 
ground state configuration. We surprisingly find, that actually the minimal 
energy is obtained when four counterions are present, hence we overcharged 
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Figure 6. Ground state configurations for two, three, four and five electrons. The corre- 
sponding geometrical figure repulsion and their typical angles are given. The electrostatic 
energy (in units of fcsTi_g/a) is given for a central charge of +2e. 



the colloid by two counterions, or by 100 %! That is, the excess counteri- 
ons gain more energy by assuming a energetically favorable configuration 
around the macroion than by escaping to infinity, the simple reason be- 
hind overcharging. In our example, the minimum is reached when four 
counterions are present. The colloid radius and the Bjerrum length enter 
as prefactors and change only the energy difference between neighboring 
states. 

The spatial correlations of the counterions are fundamental to obtain 
overcharging. Indeed, if we apply the same procedure and smear Z counte- 
rions onto the surface of the colloid of radius a, we obtain for the energy 



E = l B 



1 Z^ Zrn Z 

— . 29 

2 a a 



The minimum is reached for Z = Z m , hence no overcharging can occur. 

The important message to be learned is that, from an energetic point of 
view, a colloid always tends to be overcharged by discrete charges. Other 
important geometries like infinite rods or infinitely extended plates cannot 
be treated in such a simple fashion because they are not finite in all direc- 
tions. One needs therefore enough screening charges in the environment to 
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Figure 7. Electrostatic energy (in units of fcsTb) for ground state configurations of a 
single charged macroion as a function of the number of overcharging counterions n for 
three different bare charges Z m . The neutral case was chosen as the potential energy 
origin, and the curves were produced using the theory of Eq. (33), compare text. 



limit the range of the interactions in the infinite directions, there is a need 
for a minimal amount of salt present to allow for overcharging[102], which 
is not the case for a colloid. 

Obviously, for a large number of counterions the direct computation of 
the electrostatic energy by using the exact equation (28) becomes unfeasi- 
ble. Therefore we resort to simulations for highly charged spheres. 

5.2. ONE COLLOID 

The electrostatic energy as a function of the number of overcharging coun- 
terions n is displayed in Fig. 7. We note that the maximal (critical) accep- 
tance of n (4, 6 and 8) increases with the macroionic charge Z m (50, 90 
and 180 respectively). Furthermore for fixed n, the gain in energy is always 
increasing with Z m . Also, for a given macroionic charge, the gain in energy 
between two successive overcharged states is decreasing with n. 

In the ground state the counterions are highly ordered. Rouzina and 
Bloomfield [91] first stressed the special importance of these crystalline 
arrays for interactions of multivalent ions with DNA strands, and later 
Shklovskii [81] showed that the Wigner crystal (WC) theory can be applied 
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to determine the interactions in strongly correlated systems. In two re- 
cent short contributions [106, 107] we showed that the overcharging curves 
obtained by simulations of the ground state, like Fig. 7, can be simply 
explained by assuming that the energy e per counterion on the surface 
of a macroion scales as y/c, where c denotes the counterion concentration 
c = N/A, N is the total number of counterions on the surface and A the 
total macroion area. This can be justified by a simple argument, where 
each ion interacts in first approximation only with the oppositely charged 
background of its Wigner-Seitz (WS) cell, which can be approximated by 
a disk of radius h, yielding the same WS cell area. 

For fixed macroion area we can write the energy per counterion as 

£ W(N) = -^SVN = -aWiyfi, (30) 
VA 

where £ = IbZ 2 and the simple hole theory gives = 2^/tt 3.54[108]. 

For an infinite plane, where the counterions form an exact triangular lat- 
tice, one obtains the same functional form as in Eq. (30), but the prefactor 
a^> gets replaced by the numerical value a wc = 1.96[109]. 

Not knowing the precise value of a we can still use the simple scaling 
behavior with c to set up an equation to quantify the energy gain AE\ by 
adding the first overcharging counterion to the colloid. To keep the OCP 
neutral we imagine adding a homogeneous surface charge density of op- 
posite charge ( ~^ c6 ) to the colloid [81]. This ensures that the background 
still neutralizes the incoming overcharging counterion and we can apply Eq. 
(30). To cancel our surface charge addition we add another homogeneous 
surface charge density of opposite sign This surface charge does not in- 
teract with the now neutral OCP, but adds a self-energy term of magnitude 
\ | , so that the total energy difference for the first overcharging counterion 
reads as 

A .Ed = (N c + l)e(N c + 1) - N c e{N c ) + A (31) 
By using Eq. (30) this can be rewritten as [110] 
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Completely analogously one derives for the energy gain AE n for n over- 
charging counterions 
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Using Eq. (33), where we determined the unknown a from the simulation 
data for AE\, we obtain a curve that matches the simulation data almost 
perfectly (Fig. 7). The second term in Equation (33) also shows why the 
overcharging curves of Fig. 7 are shaped parabolically upwards for larger 
values of n. If one successively removes each of n counterions from a neutral 
colloid, one can derive in a similar fashion the ionization energy cost 



Using the measured value of a we can simply determine the maximally ob- 
tainable number n ma x of overcharging counterions by finding the stationary 
point of Eq. (33) with respect to n: 



The value of n max depends only on the number of counterions N c and a. 
For large N c Eq. (35) reduces to -j^7= y/Nc which was derived in Ref . 

[105] as the low temperature limit of a a neutral system in the presence of 
salt. What we have shown is that the overcharging in this limit has a pure 
electrostatic origin, namely it originates from the energetically favorable 
arrangement of the ions around a central charge. We also showed in Ref. 
[110] that a reaches the perfect WC value of 1.96 if the colloid radius a 
gets very large at fixed c, or when c becomes large at fixed a. 

If instead of a central charge scheme one uses discrete charge centers 
distributed randomly over the colloidal surface we find counterion struc- 
tures which are quite far away from the WC array, especially when the 
counterions are pinned to their counter charges. This depends on the in- 
teraction energy at contact, which depends of course on I and distance of 
closest approach. However, we still find overcharging, although reduced in 
value, of the form given by Eq. 33 [108, 111] 

5.2.1. Macroion- counterion interaction profile atT = OK 
The interaction profile between a completely neutralized macroion and one 
excess counterion is obtained by displacing adiabatically the excess coun- 
terion from infinity towards the macroion. From far away the counterion 
sees only a neutral object and has no measurable interaction, whereas upon 
approach to the macroion the WC hole gets created in the counterion layer, 
and we observe a distance dependant attraction towards the macroion. We 
investigated cases of Z m = 2 . . . 288. All curves can be nicely fitted with an 
exponential fit of the form 
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where AE\ is the measured value for the first overcharging counterion, and 
r is the only fit parameter. In all our results for r versus \fW c we observe a 
linear dependence for a wide range of values for N c , r oc y/N c , which again 
can be explained by applying the WC hole picture [110]. 

5.3. TWO COLLOIDS 

Now we apply what we have learned about a single colloid to two equal- 
sized, fixed charged spheres of bare charge Qa and Qb separated by a 
center-center separation R and surrounded by their neutralizing counteri- 
ons, which give concentrations ca and cb, respectively. 

All these ions making up the system are immersed in a cubic box of 
length L = 80a, and the two macroions are held fixed and disposed sym- 
metrically along the axis passing through the centers of opposite faces. This 
leads to a colloid volume fraction f m = 2- |7r(a/L) 3 ~ 8.4 x 10 -3 . For finite 
colloidal volume fraction f m and temperature, we know from the study car- 
ried out above that in the strong Coulomb coupling regime all counterions 
are located in a spherical "monolayer" in contact with the macroion. Here, 
we investigate the mechanism of strong, long range attraction stemming 
from monopole contributions; that is, one colloid is overcharged and the 
other one undercharged. 

5.3.1. Observation of metastable ionized states 

For the charge symmetrical situation we have ca = cb- When we brought 
this system to room temperature To and generated initially the counterions 
randomly inside the box we observed in some cases that one of the colloids 
remained undercharged, and the other one was overcharged, and these con- 
figurations turned out to be extremely long lived in the course of our MD 
simulations( more than 10 8 MD time steps). However it is clear that such 
a state is "metastable" because by symmetry arguments it cannot be the 
lowest energy state. The observed barrier is the result of the WC attrac- 
tion, because close to the macroion surface the energy is reduced. For very 
distant macroions the barrier height for the first overcharged state has to 
equal AE\ from Eq. 33. The barrier profile at T = can also be extremely 
well approximated by an application of Eqs. (33) and (34), plus taking 
into account the distant dependent monopole contribution [106]. This leads 
to a barrier height which scales as \fc for large separations. For smaller 
separations one has to take into account also the effect of strong mutual 
polarization of both macroions, which leads effectively to a sharing of their 
proximal counterion layer into a superlattice. This can be taken into ac- 
count by a higher effective counterion density close to the surface, leading 
to an almost linear scaling of the barrier height with c [106, 110]. 
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5.3.2. Asymmetrically charged colloids 

The most interesting phenomenon, however, appears when the two colloids 
have different counterion concentrations, here ca > cb, since then stable 
ionized states can appear. The physical reason is that a counterion can 
gain more energy by overcharging the colloid with ca then it loses by ioniz- 
ing colloid B. A straight forward application of the procedure outlined for 
the barrier calculation [107, 110] yields a simple criterion (more specifically 
a sufficient condition), valid for large macroionic separations, for the charge 
asymmetry \/Na — \/Nb to produce an ionized ground state of two unlike 
charged colloids with the same size: 



5.3.3. Finite temperature analysis 

We have also demonstrated that the ground state phenomena survive for 
finite temperatures, i.e. an ionized state can also exist at room temperature 
To. The left part of Figure 8 shows the time evolution of the electrostatic 
energy of a system Za = 180 with Zb = 30, R/a = 2.4 and a colloidal 
volume fraction of 7 • 10~ 3 , where the starting configuration is the neutral 
state {DI = 0). One clearly observes two jumps in energy, AE\ = —19.5 
and AE2 = —17.4, which corresponds each to a counterion transfer from 
colloid B to colloid A. These values are consistent with the ones obtained 
for the ground state, which are— 20.1 and —16.3 respectively. Note that this 
ionized state {DI = 2) is more stable than the neutral but is expected to be 
metastable, since it was shown previously that the most stable ground state 
corresponds to DI = 5. The other stable ionized states for higher DI are 
not accessible with reasonable computer time because of the high energy 
barrier made up of the correlational term and the monopole term which 
increase with DI. In the right part of Fig. 8 we display a typical snapshot 
of the ionized state {DI = 2) of this system at room temperature. 

Obviously, these results are not expected by the DLVO theory even in 
the asymmetric case (see e. g. [112]). Previous simulations of asymmetric 
(charge and size) spherical macroions [113] were also unable to predict such 
a phenomenon since the Coulomb coupling was weak (water, monovalent 
counterions). Note that the appearance of (meta-)stable ionized states can 
alter the effective interactions between charged colloids in solution. The 
monopole attraction will lead to attraction between like charged colloids, 
flocculation, and related phenomena. 

At this stage, we would like to stress again, that the appearance of 
a stable ionized ground state is due merely to correlation. An analogous 
consideration with smeared out counterion distributions along the lines of 
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Figure 8. Relaxation, at room temperature To = 298-/f, of an initial unstable neutral 
state towards ionized state. Plotted is the total electrostatic energy versus time (LJ 
units), for Zb = 30 and R/a = 2.4. Dashed lines lines represent the mean energy for each 
DI state. Each jump in energy corresponds to a counterion transfer from the macroion 
B to macroion A leading to an ionized state that is lower in energy than the neutral one. 
The right figure is a snapshot of the final ionized state, with net charges +4e and ~4e as 
indicated. 

Eq. (29) will again always lead to two colloids exactly neutralized by their 
counterions [114]. Our energetic arguments are quite different from the 
situation encountered at finite temperatures, because in this case even a PB 
description would lead to an asymmetric counterion distribution. However, 
in the latter case this happens due to purely entropic reasons, namely in the 
limit of high temperatures, the counterions want to be evenly distributed 
in space, leading to an effective charge asymmetry. 

Note also, that there can exist parameter regions, such as high molar 
electrolytes, where the overcharging of a single macroion is due to mainly 
entropic effects [98, 115, 116], whose exact mechanism is currently under 
investigation [11 7] . 
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